############################################################################
#Plot for marginal effects figures
#Irene Menendez, March 2016
#To correct for over-confidence, the code follows Esarey and Sumner (2015)
#and uses their package interactionTest
############################################################################

library(interactionTest)

############################################
#First: Marginal effect of nhhi2 on lgenben
############################################

##(a) Low district magnitude
#Get marginal effects and standard error from STATA's margins command
#This is from following model
#quietly xi:xtreg lgenben c.mgdp##c.dm2##c.nhhi2 xgdp unemp capgdp pop15_64 lgenbenL1, fe cluster(cntry)
#margins, dydx(mgdp) at (nhhi2=(0.004(0.07)0.38) dm2=(1))
#Vector with values of conditioning variable nhhi2
z <- c(.004, .044, .084, .124,  .164,  .204,  .244, .284,  .324, .364)

#Vector of marginal effects
m1 <- c(-.017812,  -.0119823, -.0061526,  -.0003229, .0055069, .0113366, .0171663,  .022996, .0288257, .0346554)

#Vector of standard errors
s1 <- c(.0076992, .0058442,.004254, .0033323, .0036299, .0049309, .0066721, .0085898,  .0105885, .0126299)

#Degrees of freedom
df1 <- 351

#All these are input for the function
tc1 <- fdrInteraction(m1, s1, df1, level = 0.90)

#Calculate lower bound of corrected CI and original CI
lbc1 <- m1 - tc1*s1
lbo1 <- m1 - 1.65*s1

#Calculate upper bound of corrected CI
ubc1 <- m1 + tc1*s1
ubo1 <- m1 + 1.65*s1

#par(mfrow=c(1,1))
plot(z, m1, ylim = c(-0.04, 0.06),
pch = 19, col = "black",
main = "Low district magnitude",
xlab = "Geographical Concentration", # (i.e., nhhi),
ylab = "Marginal effect of imports on benefit generosity",
cex.lab = 1.1
)

segments(z, lbc1, z, ubc1)
abline(h=0, lty = "dotted")

#Indicate 90% CI without adjustment for multiple testing by horizontal bars
segments(z - 0.003, lbo1, z + 0.003, lbo1)
segments(z - 0.003, ubo1, z + 0.003, ubo1)

##(b) High district magnitude
m2 <- c(.0478266, -.0067347, -.0612959, -.1158572, -.1704185, -.2249797, -.279541, -.3341023, -.3886636,  -.4432248)

s2 <- c(.0212475, .0048768, .025349, .0480398, .0708486, .0936894, .1165432, .1394038, .1622682 , .185135)

tc2 <- fdrInteraction(m2, s2, df1, level = 0.90)

#lower bound of corrected CI and original CI
lbc2 <- m2 - tc2*s2
lbo2 <- m2 - 1.65*s2

#upper bound of corrected CI
ubc2 <- m2 + tc2*s2
ubo2 <- m2 + 1.65*s2

#ylim = c(-0.6,0.1),

plot(z, m2, ylim = c(-0.8,0.2),
pch = 19, col = "black",
main = "High district magnitude",
xlab = "Geographical Concentration", # (i.e., nhhi),
ylab = "Marginal effect of imports on benefit generosity",
cex.lab = 1.1
)

segments(z, lbc2, z, ubc2)
abline(h=0, lty = "dotted")

#Indicate 90% CI without adjustment for multiple testing by horizontal bars
segments(z - 0.003, lbo2, z + 0.003, lbo2)
segments(z - 0.003, ubo2, z + 0.003, ubo2)

#############################################
#Second: Marginal effect of nhhi2 on lgenalmp
#############################################

##(a) Low district magnitude
z <- c(.004, .044, .084, .124,  .164,  .204,  .244, .284,  .324, .364)
m3 <- c(-.0110928, -.0018079, .0074771, .016762, .026047, .035332, .0446169,  .0539019, .0631868, .0724718)

s3 <- c(.0156402, .0132007, .0115881, .0111664,  .0120612, .0140228, .0166789, .0197513, .0230742, .0265539)

df3 <- 307

tc3 <- fdrInteraction(m3, s3, df3, level = 0.90)

#lower bound of corrected CI and original CI
lbc3 <- m3 - tc3*s3
lbo3 <- m3 - 1.65*s3

#upper bound of corrected CI
ubc3 <- m3 + tc3*s3
ubo3 <- m3 + 1.65*s3


plot(z, m3, ylim = c(-0.05,0.15),
pch = 19, col = "black",
main = "Low district magnitude",
xlab = "Geographical Concentration", # (i.e., nhhi),
ylab = "Marginal effect of imports on ALMP generosity",
cex.lab = 1.1
)

segments(z, lbc3, z, ubc3)
abline(h=0, lty = "dotted")

#Indicate 90% CI without adjustment for multiple testing by horizontal bars
segments(z - 0.003, lbo3, z + 0.003, lbo3)
segments(z - 0.003, ubo3, z + 0.003, ubo3)

##(b) High district magnitude
m4 <- c( .0750113, .0180893, -.0388327, -.0957547, -.1526768, -.2095987, -.2665208, -.3234428, -.3803647, -.4372867)

s4 <- c( .0317294,  .0129468, .0449813, .0810793, .1175552, .1541409, .1907733, .2274299, .2641005, .30078)

tc4 <- fdrInteraction(m4, s4, df3, level = 0.90)

#lower bound of corrected CI and original CI
lbc4 <- m4 - tc4*s4
lbo4 <- m4 - 1.65*s4

#upper bound of corrected CI
ubc4 <- m4 + tc4*s4
ubo4 <- m4 + 1.65*s4

plot(z, m4, ylim = c(-1.2,0.5),
pch = 19, col = "black",
main = "High district magnitude",
xlab = "Geographical Concentration", # (i.e., nhhi),
ylab = "Marginal effect of imports on ALMP generosity",
cex.lab = 1.1
)

segments(z, lbc4, z, ubc4)
abline(h=0, lty = "dotted")

#Indicate 90% CI without adjustment for multiple testing by horizontal bars
segments(z - 0.003, lbo4, z + 0.003, lbo4)
segments(z - 0.003, ubo4, z + 0.003, ubo4)








